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We investigate the Rubinstein-Duke model for polymer reptation by means of density-matrix 
renormalization group techniques both in absence and presence of a driving field. In the former case 
the renewal time t and the diffusion coefficient D are calculated for chains up to A*' = 150 reptons 
and their scaling behavior in A'' is analyzed. Both quantities scale as powers of A'': r ~ and 
D ~ with the asymptotic exponents z — 3 and a:: = 2, in agreement with the reptation theory. 

For an intermediate range of lengths, however, the data are well-fitted by some effective exponents 
whose values are quite sensitive to the dynamics of the end reptons. We find 2.7 < z < 3.3 and 
1.8 < X < 2.1 for the range of parameters considered and we suggest how to influence the end 
reptons dynamics in order to bring out such a behavior. At finite and not too small driving field, 
we observe the onset of the so-called band inversion phenomenon according to which long polymers 
migrate faster than shorter ones as opposed to the small field dynamics. For chains in the range of 
20 reptons we present detailed shapes of the reptating chain as function of the driving field and the 
end repton dynamics. 

PACS numbers: 83.10.Nn, 05.10.-a, 83.20.Fk 



I. INTRODUCTION 

The idea of reptation was introduced about thirty 
years ago by de Gennes j^] in order to explain some dy- 
namical properties of polymer melts of high molecular 
weight. The dynamics of such systems is strongly influ- 
enced by entanglement effects between the long polymer 
chains. The basic idea of reptation is that each polymer 
is constrained to move within a topological tube due to 
the presence of the confining surrounding polymers [^,^ . 
Within this tube the polymer performs a snake-like mo- 
tion and advances in the melt through the diffusion of 
stored length along its own contour. One focuses thus 
on the motion of a single test chain, while the rest of 
the environment is considered frozen, i.e. as formed by a 
network of fixed obstacles. The reptation theory predicts 
that the viscosity /i and longest relaxation time r (known 
as renewal time) scale as /i ~ r ~ N'^, where N is the 
length of the chains, while the diffusion constant scales 
a.s D ^ l/N"^. These results are not far from the exper- 
imental findings. Measurements of viscosity of concen- 
trated polymer solutions and melts of different chemical 
composition and nature are all consistent with a scal- 
ing ^ ~ N^-'^ [Q, while for the diffusion constant both 
D ~ 1/A^2 H and I? - l/iV2-4 have been reported. 
This discrepancy triggered a substantial effort in order 
to reconcile theory and experiments. 

Another physical situation where reptation occurs is 
in gel electrophoresis, where charged polymers diffuse 
through the pores of a gel under the influence of a driving 
electric fleld . The gel particles form a frozen network 
of obstacles in which the polymer moves through the dif- 
fusion of stored length. Electrophoresis has important 
practical applications, for instance in DNA sequencing. 



since it is a technique which allows to separate polymers 
according to their length. 

The previous examples demonstrate how reptation is 
an important mechanism for polymer dynamics in dif- 
ferent physical situations. A prominent role in under- 
standing the subtle details of the dynamics was played 
by models deflned on the lattice, which besides describ- 
ing correctly the process at the microscopic level, offer 
important computational advantages. The aim of this 
paper is to investigate in detail one of such lattice mod- 
els, which was originally introduced by Rubinstein [Tl| ] 
and later extended by Duke in order to include the 
effect of an external driving field. The Rubinstein - Duke 
(RD) model has been studied in the past using several 
techniques; there are a limited numbers of exact results 
available [^-15| while a rich literature on simulation re- 
sults on the model exists. The latter have been mostly 
obtained by Monte Carlo (MC) simulations 
MC simulations are tedious for long chains since the re- 
newal time scales as N^, thus it is difficult to obtain a 
small statistical error for large N . 

In this paper we study the RD model by means of 
Density Matrix Renormalization Group (DMRG) [p^ , 
a technique which has been quite successful in recent 
years, in particular in application to condensed matter 
problems as quantum spin chains and low-dimensional 
strongly correlated systems Using the formal sim- 

ilarity between the Master equation for reptation and 
the Schrodinger equation, DMRG also allows to calcu- 
late accurately stationary state properties of long reptat- 
ing chains. 

Some of the results reported here, in particular con- 
cerning the scaling of the renewal time r, have been pre- 
sented before pO]. Here we will give full account of the 
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details of the calculations and present a series of new re- 
sults concerning reptation in the presence of an electric 
field. One of the main conclusions of our investigation is 
that the exponents describing the scaling of r and D in 
the intermediate length region appear to be rather sen- 
sitive on the structure of the end repton. Wc find that, 
by influencing the end repton dynamics, one has a regime 
where the effective exponents are considerably lower than 
the standardly found experimental values. Experiments 
are suggested, involving polymer architectures not inves- 
tigated so far, and which ought to confirm our predic- 
tions. 

This paper is organized as follows: in Sec. || we intro- 
duce the RD model, while in Sec. HI we briefly outline 
the basic ideas of DMRG. Sec. IV is dedicated to the 



properties of reptation in absence of an external field, in 
particular to the scaling properties of r and D. Sec. 
collects a series of results of reptation in a field, while 
Sec. ^ concludes our paper. 



II. THE RUBINSTEIN - DUKE MODEL 



In the RD model the polymer is divided into N units, 
called reptons, which are placed at the sites of a d- 
dimensional hypercubic lattice (see Fig. |l|). The number 
of reptons that each site can accommodate is unlimited 
and self-avoidance effects are neglected. Each configura- 
tion is projected onto an axis along the (body) diagonal 
of the unit cell and it is identified by the relative coor- 
dinates t/i = Zi+i — Zi of neighboring reptons along the 
chain [zi indicates the projected coordinate of the i-th 
repton). The relative coordinates can take three values 
Ui = —1,0, 1 and there are thus in total 3^^^ different 
configurations for a chain with N reptons. 




FIG. 1. Example of a configuration in the RD model in 
d — 2 dimensions. In terms of the relative coordinates 
projected along the field direction the configuration reads 
j/ = {1,0, 1,1,— 1,1,0}. A non-zero field (e) biases the motion 
of the reptons, which occurs with rates B — exp{e/2) {B~^) 
for moves in the direction (opposite to) the applied field. The 
arrows indicate the possible moves. Notice that an end repton 
can stretch to d lattice positions forward and backward in the 
field (as repton 8 in the example) . 



When two or more reptons accumulate at the same lat- 
tice site they form a part of stored length, which can then 
diffuse along the chain. In terms of relative coordinates 
a segment of stored length corresponds to yi = 0, there- 
fore allowed moves are interchanges of O's and I's, i.e. 
0,±1 ^ ±1,0. On the contrary, the end reptons of the 
chain can stretch (0 ±1) or retract to the site occupied 
by the neighboring repton (±1 —^ 0). The dynamics of 
the chain is fully specified once the rates for the moves are 
given. A field e along the projection axes is introduced 
so that moves forward and backward in the field occur 
with rates B = exp(e/2) and — exp(— e/2). In the 
following we will be interested in both cases £ = and 
e > 0. Notice that when an end repton moves towards an 
empty lattice site (tube renewal process) it has in total d 
different possibilities of doing so forward and backward 
m the field (see Fig. 0), therefore the associated rates are 
dB and dB~^ respectively. On the contrary the end rep- 
ton can retract by moving to the (unique) site occupied 
by its neighbor with rate B or B^^. Summarizing the 
possible moves are: (a) stored length diffusion for inner 
segments, i.e. (0, ±1 ^ ±1, 0) with rates B and B~^, (b) 
contractions for external reptons (±1^0) with rates B 
and B"^ and (c) stretches for external reptons (0 ±1), 
with rates dB and dB~^. Notice that the parameter 
d enters in the model only at the stretching rates (c). 
Rather than linking d to the lattice coordination number 
we interpret it as the ratio between stretching rates and 
rates associated to moves of inner reptons. This allows 
to choose any positive values for d. 





FIG. 2. (a) In thick line we indicate a test chain moving 
in an environment of other chains. By modifying its ends, for 
instance, by attaching molecules of large size one can lower 
the parameter d, as in this way one expects that stretchings 
of the chain out of the tube will occur to lower rate, (b) The 
lowering of d would also occur for reptating polymers with 
short branching ends. 



As we will see, a variation of d has some important ef- 
fects on the corrections to the asymptotic behavior. This 
regime may turn out to be relevant for the experiments. 
Moreover it is conceivable that the parameter d could be 
somewhat modified in an experiment. Figure ^ schemat- 
ically illustrates a possibility of doing so, namely by at- 
taching to the chain ends some large molecules so that 
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for each chain tube renewal moves would be impeded. 
The same effects would be possible if the linear chain is 
modified such as to have branches close to the endpoints 
(Fig. §(b)), or an end part stiffer than the rest of the 
chain (as it could be realized in block copolymers). In 
both cases we expect that stretches out of the confining 
tube would be suppressed while chain retractions would 
not be much impeded. 

Once the rates for elementary processes are given, the 
stationary properties of the system can be found from 
the solution of the Master equation 



dPiy,t) 
dt 



(1) 



in the limit t ~* oo. Here P{y, t) indicates the probability 
of finding the polymer in a configuration y at time t and 
the matrix H contains the transition rates per unit of 
time between the different configurations of the chains, 
as given in the rules discussed above. 



III. DENSITY MATRIX RENORMALIZATION 

DMRG was introduced in 1992 by S. White as an 
efficient algorithm to deal with a quantum hamiltonian 
for one dimensional systems. It is an iterative basis- 
truncation method, which allows to approximate eigen- 
values and eigenstates, using an optimal basis of size m. 
It is not restricted to quantum system; it has also been 
successfully applied to a series of problems, ranging from 
two dimensional classical systems to stochastic pro- 
cesses 1^ . The basic common feature of these problems 
is the formal analogy with the Schrodinger equation for 
a one dimensional many-body system. DMRG can also 
be applied to the Master Equation (|^) for the reptat- 
ing chain, albeit with some limitations due to the non- 
hermitian character of the matrix H. In this specific case 
we start from small chains, for which H can be diagonal- 
ized completely, and construct effective matrices repre- 
senting H for longer chains through truncation of the 
configurational space followed by an enlargement of the 
chain. This is done through the construction of a reduced 
density matrix whose eigenstates provide the optimal ba- 
sis set, as can be shown rigorously (see Refs. JT8| , p^ for 
details). The size m of the basis remains fixed in this 
process. By enlarging m one checks the convergence of 
the procedure to the desired accuracy. Hence m is the 
main control parameter of the method. In the present 
case we found that m = 27 is sufficient for small driving 
fields and we kept up to m = 81 state for stronger fields. 

H is only hermitian for zero driving field. So for a fi- 
nite driving field one needs to apply the non-hermitian 
variant of the standard DMRG algorithm 1 22 1 . For a non- 
hermitian matrix H one has to distinguish between the 
right and left eigenvector belonging to the same eigen- 
value. Since iJ is a stochastic matrix the lowest eigen- 
value equals and the corresponding left eigenvector is 



trivial (see next section). The right eigenvector gives 
the stationary probability distribution. The next to low- 
est eigenvalue, the gap l/r, yields the slowest relaxation 
time T of the decay towards the stationary state. Gener- 
ally the DMRG method works best when the eigenvalues 
are well separated. For long chains and stronger driving 
fields the spectrum of H gets an accumulation of eigen- 
values near the zero eigenvalue of the stationary state. 
This hampers the convergence of the method seriously 
and enlarging the basis m becomes of little help, while 
standardly this improves the accuracy substantially. In 
order to construct the reduced density matrix from the 
lowest eigenstates one needs to diagonalize the effective 
matrices H at each DMRG step, we used the so-called 
Arnoldi method which is known to be particularly stable 
for non-hermitian problems [E3|. 



IV. ZERO FIELD PROPERTIES 

In this section we present a series of DMRG results on 
the scaling behavior as function of the polymer length 
N of the so-called tube renewal time t{N), i.e. of the 
characteristic time for reptation, and of the diffusion co- 
efficient D{N). 



A. Tube renewal time 

The renewal time, r, is the typical time of the repta- 
tion process, i.e. the time necessary to loose memory of 
any initial configuration through reptation dynamics, r 
is given by the inverse of the smallest gap of the matrix 
H, which can be seen as follows: Starting from an initial 
configuration one has, asymptotically for sufficiently long 
times {t — !■ oo): 



\P{t)) = I0o) 



(2) 



where |0o) and |(/>i) are eigenstates of the matrix H de- 
fined in Eq. (|l|) with eigenvalues Eq = and Ei = 1/t, 
respectively. The state |0o) is the stationary state of the 
process, while the gap of H corresponds to the inverse 
relaxation time. Notice that, as the matrix H is sym- 
metric at zero field, Ei is always a real number, while 
for non-Hcrmitian matrices the eigenvalues may get an 
imaginary part, which causes a relaxation to equilibrium 
through damped oscillations. 

An important point of which we took advantage in 
the calculation is the fact that the ground state |(/)o) of 
the matrix H is known exactly in absence of any driv- 
ing fields. Since H is stochastic and symmetric (i.e. 

Elyyi = '^yi Hy'y =0), it follows immediately that 
the stationary state is of the form 0o(y) — c, with c a 
constant. Instead of applying the DMRG technique di- 
rectly to the matrix H we applied it to the matrix H' 
defined by: 
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H' = H 



(3) 



Now if we choose A > _Ei the lowest eigenvalue of H' is 
equal to Ei , therefore the problem of calculating the gap 
of H is reduced to the calculation of the ground state of 
a new (non-stochastic) matrix H' . This approach is con- 
siderably more advantageous in term of CPU time and 
memory required for the program (for more details see 
Ref. g)) 

Figure g shows a plot of InT(A^) vs. IniV for the re- 
newal time as calculated from DMRG methods for vari- 
ous values of the stretching rate d and for lengths up to 
N = 150. The data have been shifted along the abscissae 
axis by an arbitrary constant. We fitted the data by a 
linear interpolation, as it is mostly done in Monte Carlo 
simulations and in experiments. The resulting slopes pro- 
vide estimates of the exponent z, which we find to be 
sensitive to the stretching rate d. For d sufhciently large 
(d > 2) we find z « 3.3, in agreement with experiments 
Q and previous Monte Carlo simulation results pl] , p5[ . 
The value of z decreases when d is decreased. For d suf- 
ficiently small {d w 0.1) we find z w 2.7 < 3, which is a 
regime that has not yet been observed in experiments. 
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d = 0.1, z = 2.70 
d = 0.5, z = 3.00 
d= 1.0, z = 3.21 
d = 3.0, z = 3.29 
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FIG. 3. Plot of \iit{N) vs. InN for various d. We report 
the values of the slope of the data obtained from a linear 
interpolation. 

To shed some more light into these results we consid- 
ered the effective exponent^ i.e. the following quantity: 



ZN 



lnT(A^ + l)-lnT(iV-l) 
ln(7V + 1) - ln(7V - 1) 



(4) 



which is the discrete derivative of the data in the log-log 
scale plot. Such quantity probes the local slope (at a 
given length TV) of the data of Fig. ^ and provides a bet- 
ter estimate of the finite N corrections to the asymptotic 
behavior. The calculation of zjv requires very accurate 
data and cannot be easily performed by Monte Carlo sim- 
ulation results which are typically affected by numerical 



uncertainties, as these are amplified when taking numer- 
ical derivatives. We stress that already in the log-log 
plot of Fig. H the deviation of the data from linearity 
is noticeable, therefore the values of z given above are 
just average values and not to be expected as asymptotic 
ones. 

Figure || shows a plot of zat for the data given in Fig. 

plotted as function of 1/y/N. In the thermodynamic 
limit ^ cx), zjv is seen to converge towards z — i, 
in agreement with de Gennes' theory [|l|. Corrections to 
this asymptotic limit yield zjv > 3 when N is sufficiently 
large, with deviations towards zat < 3 for small d and 
not too large N . 
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FIG. 4. Solid lines: Effective exponent zn for various d. 
Dashed line: Corrections due to CLF as predicted by Doi 
(Eqs. ^ and ^) for the d = 3 curve. Dot-dashed line: In- 
clusion of higher order corrections as given by Eq. ^ for the 
curve d = 1. 

There exist some theoretical predictions for the form of 
the corrections to the asymptotic behavior of the renewal 
time . These are based on contour length fluctuations 
(CLF), i.e. on the idea that the length of the tube fluc- 
tuates in time and that this would help accelerating the 
renewal process. Such fluctuations were not taken into 
account in the original work of de Gennes, who assumed 
the tube to have a fixed length. According to CLF the- 
ory, thus, T scales as [p6l 



t{N) ~ iV^ 



1 - 




(5) 



with A'^o a typical length. Substituting this equation into 
Eq. (0) one obtains: 



JCLF) 



(6) 



which is a monotonically increasing function of 1/ \/N. 
For all lengths in the physical range N > Nq one has 
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z^^^^ > 3. Eq. @ is plotted as a dashed line in Fig. 
0, where we have chosen Nq = 2.3 in order to obtain 
the best fit of the DMRG data with d = 3. This value 
is consistent with that predicted by the CLF theory 
Our results further suggest that the value of A^o increases 
when d is decreased. 

While for sufficiently long chains the data apparently 
merge to the CLF theory given by Eq.(^ higher order 
corrections appear to be of opposite sign. When d is 
lowered the latter become particularly strong so that the 
effective exponent for a certain range of lengths is even 
smaller than 3. To our knowledge such an effect has 
not yet been observed in experiments. Presumably stan- 
dard polymer mixtures will correspond to d ^ 1, which 
is a regime where an exponent z « 3.4 is observed [Q. 
It is conceivable that mixtures with modified architec- 
ture, as those illustrated in Fig. ||(b), i.e. long polymers 
with short branching ends, would correspond to curves 
at much lower d. Experimental results for such systems 
would be of high interest in order to check the predictions 
of the RD model in the low d regime. 

In order to gain some more insight in the precise form 
of t{N) we considered higher order terms which lead to 
an expression of the type: 



r(iV) - 




A- 



N 



(7) 



Recently, Milner and McLeish [g^ formulated a more 
complete theory beyond that of CLF by Doi, using ideas 
from the theory of stress relaxation for star polymers [ p8| , 
which to lowest orders in 1/Vn yields an expression of 
the type given above. The effective exponent now reads: 



ZN 



No l-il + A)^/m/N 



(l-y/m/N^ +ANo/N 



(8) 



Notice that the CLF expression (Eq. ^) diverges for 
TV Nq, while the previous formula this divergence 
does not occur. For N not too large the numerator in 
Eq. (^) may change sign, reproducing features found in 
the DMRG calculations, i.e. an effective exponent z < 3. 
The dot-dashed line of Fig. ^represents a fit for the d = 1 
case using Eq.(|8|); we find that the choice Nq — 3.6 and 
A = 0.44 fits very well the numerical data for N > 25, 
while deviations for shorter chains are clearly visible. We 
stress that Eq. (|^) fits quite well the renewal time data 
in the cage model which is another lattice model 
of reptation dynamics [|T]. Finally, very recently the ef- 
fect of constraints release (CR), introduced in the RD 
model in a self-consistent manner, has been considered 
pof . Apart from small quantitative shifts in the effective 
exponents, the is unchanged. 



B. Diffusion coefficient 



We consider now a similar scaling analysis of the diffu- 
sion coefficient D. To compute the diffusion constant as 
function of the chain length we used the Nernst-Einstein 
relation ITl: 



D = lim . 

e^o Ne 



(9) 



where v is the drift velocity of the polymer with N rep- 
tons and subject to an external field e. In order to esti- 
mate D{N) we considered small fields (down to e ~ 10"^) 
and calculated the velocity Vi of the j-th repton in the 
stationary state using the formulas given in Section |VB| . 
As a check for the accuracy of the procedure we verified 
that the drift velocity Vi is independent on i, i.e. the 
position of the repton along the chain in which it is com- 
puted. 

TABLE I. Comparison between diffusion coefficients from 
exact diagonalization methods, as given in Ref. and as 
obtained from the calculation of r(e) — v/Ne from the DMRG 
algorithm with m = 18 for decreasing e. 



N 


r(e = 10-") 


r{e = 10-3) 


r{e = 10-*) 


DN^ (Ref. [p7|) 


5 


0.864873 


0.864907 


0.864908 


0.864908 


7 


0.769024 


0.769057 


0.769057 


0.769057 


9 


0.703975 


0.703951 


0.703951 


0.703951 


11 


0.657701 


0.657550 


0.657549 


0.657549 


13 


0.623372 


0.623010 


0.623006 


0.623006 


15 


0.596975 


0.596304 


0.596297 


0.596297 


17 


0.576097 


0.575007 


0.574996 


0.574996 


19 


0.559224 


0.557595 


0.557580 


0.557579 
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FIG. 5. Log-log plot of the diflusion constant as function 
of the chain length for various values of the stretching rate d. 
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In practice we estimated the ratio r = vN/e for de- 
creasing e until convergence was reached (typically we 
considered e ~ 10^** — lO^'^). Table | shows the numeri- 
cal values of r{e) for various e obtained from DMRG with 
m = 18 states kept, for short chains, so that they can be 
compared with the exact diagonalization data shown in 
the last column. The results for e ~ 10^'* are in excel- 
lent agreement with exact diagonalization data for DN^ 
(from Ref. reported in the last column, which indi- 
cates that the procedure used to extrapolate the diffusion 
coefficient from the Nernst-Einstein relation (Eq. ||) is 
correct and reliable. 

According to reptation theory the diffusion coefficient 
for large should scale as D{N) - |]. This re- 

sults has been also derived rigorously for the RD model 
where the coefficient of the leading term is also known 

i2|,ra: 



DiN) 



1 



(2rf+ 1)7V2 



(10) 



Early experimental results on diffusion coefficients for 
polymer melts were consistent with a power — 2 while 
more recent results suggest that such exponent would be 
significantly higher 1^ . In fact the original measurements 
consistent with a scaling raised quite some prob- 

lems in the past, which was pointing to an inconsistence 
between the scaling of D and t due to the following ar- 
gument. For a reptating polymer, D and r the typical 
radius of gyration should scale as: 



(11) 



Now as the polymer obeys gaussian statistics Rg ^ vN, 
which implies that, if r ~ N^^'^ then necessarily D 
iV~2'^, for the same range of lengths. This argument is 
however not quite correct in this form since Rg ^ ^/N 
only asymptotically, thus finite size corrections may af- 
fect D and r differently, as indeed happens in the RD 
model. 

Figure || shows a plot of In _D as a function of In TV for 
various values of the stretching rate d. The best fitting 
parameters x for the scaling D ~ are given. As for 

the diffusion constant the exponent passes the presumed 
asymptotic value of 2 for d = 1 and d = 3, while x < 2 
for d = 0.5 and d — 0.1. Other numerical investigations 
of the RD and related model yielded x « 2.0 ||ll[], x « 2.5 
p5[ , and x ~ 2.0 Again, it is best to analyze the 

effective exponent 



Xn 



lnD{N +1) -liiD{N -1) 
\n{N + 1) - \n{N - 1) 



(12) 



which is shown in Fig. xn shows a similar behav- 
ior as the renewal exponent. However, comparing the 
same values of N and d in Figs. ^ and ^ one notices 
that finite A^ corrections are weaker for the diffusion co- 
efficient. For instance, for d = 3 and A^"^/^ — 0.3 one 
has zn « 3.8, while xn ~ 2.3. As mentioned above. 



recent experimental investigation of diffusion coefficient 
for polymer melts and solutions yielded some different 
results concerning its scaling behavior. Earl'v 



measure- 



|5[, while in 
§. Very 



ments of polymer melts yielded D ~ A^^^ 
concentrated solutions typically D ~ A^^^-^ 
recently, however, on hydrogenated polybutadiene con- 
centrated solutions and melts it was found D ^ A^"^ '* 
for both cases A reanalysis of previous experimental 
results on several different polymers lead to the conclu- 
sion D ~ A^"^ '^ Thus the issue experimentally has 
not yet fully settled. A theoretical analysis of the con- 
tour length fluctuations on the diffusion coefficient was 
recently performed p3] leading to the expression: 



Dclf{N)^N- 




(13) 



In Fig. ^ we plotted the corresponding effective expo- 
nent, fitted to the d = 3 results. Again, the numerical 
results seem to approach the CLF formula only for very 
long chains, where the effective exponent is already quite 
close to the asymptotic value. Our results indicate some 
variation of the effective exponent as function of the pa- 
rameter d. It would be therefore interesting to investigate 
the scaling of D for the polymer with short branching 
ends, as those illustrated in Fig. |2|(b) in order to test the 
validity of our predictions at smaller d. 
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FIG. 6. Effective exponent xn as defined by Eq. for 
various d. The dotted line is the theoretical prediction from 
Eq. 

Finally we mention that the further order correction 
terms for the scaling of D have raised recently some 
debate about the nature of the expansion [|l6|,^. The 
DMRG results, which are accurate enough to investigate 
the higher orders, clearly show that these results can be 
very well represented by a series in 1/Vn pO|. 
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V. REPTATION IN THE PRESENCE OF AN 
ELECTRIC FIELD 

For finite values of e severe limitations appear prevent- 
ing us to extend the calculations to the long chains which 
we were able analyze in the limit e ^ 0. These limita- 
tions are intrinsic to the problem and the Arnoldi rou- 
tine for finding the lowest eigenvalue of the matrix H. 
For finite e and longer chains a massive accumulation of 
small eigenvalues near the ground state eigenvalue starts 
to emerge. The problem is also present in the straight 
application of the routine to small chains, for which an 
exact diagonalization of the matrix H can be performed. 
For > 11 and e > 1 the method for finding the lowest 
eigenvalue is no longer convergent. The DMRG method, 
described above, yields convergent results up to the chain 
TV = 15. For smaller fields, e < 1, which is however the 
most interesting region from the physical viewpoint, the 
calculation can be extended to somewhat longer chains 
i.e. iV « 30 — 40. Neither extension of the truncated 
basis set, nor the inclusion of more target states is a 
remedy for the lack of convergence. Despite the intrinsic 
difficulty of treating long chains, substantial information 
can be obtained by the analysis of the drift velocity for 
intermediate chain lengths as function of e. 



A. Drift velocity vs. field 

In Fig. 1^ we have plotted the drift velocity as func- 
tion of e for chains up to = 15 and stretching ratio 
d = 1. The general behavior is a rise turning over into 
an exponential decay. For small e the rise is linear, in 
agreement with the results presented for the zero field 
diffusion coefficient. In fact we have calculated the diffu- 
sion coefficient as the d erivat ive of the drift velocity with 
respect to e (see Sec. (IVB)). The exponential decay is 
in accordance with the behavior predicted by Kolomeisky 
p5[ . For large fields the probability distribution narrows 
down to a single U shaped configuration with equally 
long arms at both sides of the center of the chain. For 
odd chains Kolomeisky shows that, in the limit e — s- oo 



exp 



N 



(14) 



For the chains up to = 7 we could confirm this be- 
havior (see inset of Fig. Q), but for the longer chains one 
would have to go to larger values of e than can be handled 
with the Arnoldi method to see the asymptotic behavior. 
Generally the curves are a compromise between two ten- 
dencies: the expression for the drift grows with the field 
but the probability on a configuration shifts towards con- 
figurations with the least velocity. With increasing e ini- 
tially the first effect dominates but gradually the second 
effect overrules the first. 
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FIG. 7. The velocities as a function of an electric field for 
chains of different lengths for d = 1. Inset: Inu vs. e for 
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shown as a dashed line. 



We have also investigated the behavior of the velocity 
as function of the stretching ratio d. In Fig. || we have 
plotted the various curves for chain length = 7 for d 
ranging from 0.05 to 3.25. We see that the overall veloc- 
ity becomes small for small d which is simply a sign of 
slowing down due to the slowed down motion of the end 
reptons. For the higher values of d we see again a decrease 
in the drift velocity, which is explained by stretching of 
the chain as we will see. 
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FIG. 8. The velocities as a function of an electric field for 
difi'erent d. The number of reptons is fixed at A'' = 7. 

A closer inspection of the curves v vs. e in the interval 
< e < 1 for the longer chains shows a definite deviation 
of the linear rise to larger values of the drift velocity. This 
interesting intermediate regime was also investigated by 
Barkema et al. U& by means of Monte Carlo simulations. 
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They propose to describe the drift velocity in this region 
by the phenomenological crossover expression 



v{e,N) 



(2d+l)N 



[l + A(e7V)2)] = 



(15) 



which seems to fit quite well the Monte Carlo data [ p6[ . 
Thus in the regime where e is small but the combination 
eN of order unity, the drift velocity becomes independent 
of N (band collapse) and quadratically dependent on e. 
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FIG. 9. Plot of the drift velocity v as function of the in- 
verse chain length for various fixed values of the electric 
field e and d = 1. Inset: Blow up of v vs. 1/A'' for e = 0.2; a 
minimum in the velocity can be clearly seen. 

In Figure ^ we plotted the velocity as function of the 
inverse chain length 1/iV for various values of the elec- 
tric field. The deviations from the linear regime, where 
V ^ s/N, are clearly observable as for sufficiently long 
chains and not too small fields the drift velocity ap- 
proaches a non- vanishing value. A closer inspection on 
the curves reveals that the limiting constant velocity in 
the asymptotic regime iV — s- oo is reached through a min- 
imum in the velocity at a given polymer length A^min- In 
our calculations this occurs for any values of the fields 
in the range e > 0.10. Presumably the minimum is 
shifted to lengths N much longer than those which can 
be reached in the present calculation. The existence of a 
minimum velocity at finite length can be clearly seen in 
the inset of Fig. || which plots v vs. 1/N for e — 0.2, but 
it is present also in other curves. According to our er- 
ror estimates typical uncertainties in the velocities are at 
the sixth decimal place, thus error bars are much smaller 
than all symbol sizes in the figure and inset. 

The possible existence of a velocity minimum has given 
rise to quite some discussions in the past and it is a phe- 
nomenon referred to as band inversion (for a review see 
Ref. [|lO|). In the band inversion region long polymers 
migrate faster along the field direction than short ones, 
a physical situation that looks at first sight somewhat 
counterintuitive. Band inversion was predicted first in 



the context of the biased reptation model (BRM) [ p7[ , 
which provides a simplified picture of the physics of poly- 
mer reptation in a field, in which fluctuations effects are 
neglected. More refined analytical calculations in mod- 
els where fluctuation effects are taken into account still 
predict the presence of a velocity minimum ||3^ ] , thus the 
original result of the BRM was confirmed. 

It should be pointed out that the Eq. ([l^) does not 
yield any minimum in the velocity as dv/dN < for 
all A^, thus it cannot be strictly correct. Although our 
results are limited to the onset of the minimum, we sus- 
pect that the minimum is not a very pronounced one so 
that for longer chains the curves v vs. 1/N would appear 
rather flat. A shallow minimum could have been thus 
missed in the Monte Carlo simulations of Ref. . 

Indeed, the data presented in Ref. |^ yield, for in- 
stance, V « 0.010 for e = 0.2 and A^ = 100, thus the 
limiting velocity appears to be approached very slowly 
from below (see inset of Fig. Other previous simu- 
lations by Duke ||l^ of the RD model yielded instead a 
quite clear velocity minimum, but it is difficult to judge 
the quality of his data as no error bars are reported. 
Probably the inversion effect is somewhat weaker than 
reported in |lj]. We investigated further the effect of 
d and found that for higher d the minimum appears to 
be somewhat more pronounced, which may explain the 
results of Ref. ||l^, where d = 6 was taken (the lattice 
coordination number for an fee lattice). 

Unfortunately, present DMRG computations are lim- 
ited to the onset of the inversion phenomenon. The prob- 
lem with the DMRG approach is that it builds up an 
optimal basis for the stationary state using information 
of stationary states for shorter chains. Around the band 
inversion point there is a kind of phase transition from 
short unoriented polymers to long oriented ones (see for 
instance Ref. (lOj). This change implies that the opti- 
mal basis for short chains may be no longer a good one 
when longer chains are considered. We alleviated some- 
what this problem using the DMRG algorithm at fixed 
A^e, which allows indeed to study slightly longer chains, 
yet not enough to go deep into the band collapse regime, 
but sufficient to bring in evidence the velocity minimum. 

Originally, it was thought that the drift velocity should 
be described by a scaling function in terms of the com- 
bination A^e^ 1 39 4^. It is nowadays believed that the 
correct scaling form for the velocity should be given by 
an expression [pd|Jl6| 



v[e,N)^- g{Ne) 



(16) 



with g{x) go > for a; ^ and g{x) ~ x for a; ^ oo 
in order to match the known behavior of the velocity at 
large and small fields. The condition to have band inver- 
sion is dv/dN — for some A^ at fixed fields, which is 
equivalent to the requirement g{x) = xg'{x), for a non- 
zero value of a; = A^e. Notice that choosing the most gen- 
eral scaling form v{e,N) ~ g{Ne°') and from the re- 
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quirement of dv/dN = one still obtains g{x) = xg'{x), 
with X = Ne". 

The scaling behavior of the minimum as function of 
the field may be used to test the velocity formula. We 
calculated the polymer length N^i^ at which the mini- 
mum of the velocity occurs as function of the applied field 
e. If Eq. (|l^ ) is correct we expect e ~ 1/Nynin- Figure 
|To| shows a plot of Ine vs. In A'min for d = 1. The data 
show some curvature due to corrections to scaling and 
approach, for the longest chains analyzed, the asymptote 
£ ~ N~°', with a = 1.4, as illustrated in the figure. As 
iVinin is increased one observes a systematic decrease of 
the local slope of the data suggesting that the asymptotic 
exponent a < 1.4. In an attempt to reach the asymptotic 

time we extrapolated the local slopes of the data in Fig. 
in the limit A^min — *■ oo, which yield an extrapolated 
value a ~ 1.1, not far from the prediction of Eq.(16), 
but not fully consistent with it. To prove Eq.([l6|) more 
convincingly one would need to investigate longer chains. 
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FIG. 10. Plot of Ine vs. IniVmin for d = 1. 



B. Correlations and profiles 

A more detailed insight in the shapes of the config- 
urations is obtained by plotting averages of the local 
variables i/i. The DMRG procedure leads naturally to 
the determination of the probabilities for two consecu- 
tive segments 



(17) 



The probabilities on a single segment follow from these 
values 

y' y' 
which in turn arc normalized 



(19) 



The 9 possible values oipi{y, y') are restricted by these 
conditions. One finds another set of relations between 
these quantities by summing the Master Equation over 
all but one segment value. Exclusion of an internal seg- 
ment yi from the summation yields for yi ^ y = ±1 



B«p,_i(0,y)-B-V<-i(y,0) = 
Byp,iO,y)-B-yp,iy,0) = viy). 



(20) 



Note that v{y) is independent of the index i of the seg- 
ment under consideration. The two relations ( |20| ) are an 
expression of the fact the average velocity of the reptons 
in the field direction and the curvilinear velocity are con- 
stant along the chain. Taking the segment value yj — 
one finds 



z;(l)-t;(-l) = 2z;(l), 



(21) 



with V the drift velocity. Excluding the end segments 
from the summation yields the 2 equations (y = ±1) 



viy) = Bypi{y)-dB-ypi{0) 
viy) = dBypNiO) ' B-ypN{y) 



(22) 



One observes that the 3 probabilities on the end seg- 
ments are fixed by the normalization and the 2 equations 
( p^ ) . In general these relations show how delicate the de- 
velopment of the correlations is. In the fieldless case the 
probability of a configuration factorizes in a product over 
probabilities of segments. So for d = 1 the probability 
on any segment becomes equal to 1/3. Then of course 
the velocities vanish. Considering the terms linear in the 
field (or in _B — B~^) one sees that the drift velocity of 
the first repton, given by 

V = Bpi(l) - V(-l) + d{B - S-i)pi(O), (23) 

requires a delicate compensation in the linear deviations 
of the probabilities pi{y) in order to give a value which 
vanishes as 1/A^ for long chains. 

Rather than giving the values oi pi{y) we plot the av- 
erages 



{y^) = P^{l) - P^{-l) 



(24) 



and 

{ytVi+i 



:p,(i,i)+p,(-i,-i)-p,(i,-i)-p,(i-,i). 

(25) 
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FIG. 11. The profiles for e = 0.001 plotted as function 
of the reduced distances /jv(j) — {2i — N)/{N — 2) and 
fcjv(i) = {2i -N+1)/(N - 3) (with i = 1, 2, . . . TV - 1). The 
number of reptons is fixed at A'' = 15. 




FIG. 12. The shapes of the chain for e — 0.001 and various 
values of d. The number of reptons is fixed at A'' = 15. 



A typical plot for the small e regime is given in Fig. 
[lT| . The global symmetry due to the interchange of head 
and tail makes the average (|^) antiysymmetric with re- 
spect to the middle and the average ( p5|) symmetric. One 
observes that the features increase with the mobility d of 
the end reptons. The averages give information about 
the average shape of the chain. In Fig. |l^ we translate 
the averages ( p4| ) into average spatial configurations by 
integrating (summing) the segment values to positions 
with respect to the middle repton. Clearly the develop- 
ment of the U shape is visible with increasing d. The 
effect is larger at the ends than in the middle. 
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FIG. 13. As in Fig. |ll| for d = 1 and varying e. 
number of reptons is fixed at A = 15. 



The 



As an example of the behavior in the intermediate 
regime (of the velocity profiles Fig. 3) we have plotted in 
Fig. ^ the situation for = 15, d = 1 and various val- 
ues of e. For the larger values of e the average (j/i) does 
not change very much but the plateau of the correlations 
{yiUi^i) in the middle keeps rising. We expect that for 
longer chains a large region in the middle develops where 
the shape varies weakly with e and where the correlations 
between consecutive segments increase. Thus the chain 
obtains longer stretches which are oriented in the field, 
either up or down, but which largely compensate, such 
that the overall shape in the middle remains more or less 
the same. This is in agreement with the speculations of 
Barkema et al. |l^ on the chain as a stretched sequence 
of more or less isotropic " blobs" . 
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FIG. 14. The profiles for d = 1. The number of reptons is 
fixed at A = 9. 
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Finally we show in Fig. |lj the situation for strong e 
on a chain of iV = 9 reptons. For the strongest values of 
e it is almost exclusively in the U shaped configuration. 
Note also that the correlations approach 1 except for the 
middle pair of segments which arc on different branches 
of the U. So ultimately the value of the correlation in 
the middle will approach -1. 

Thus we see that the reptating chain develops a rich 
and delicate pattern of shapes and correlations which are 
not so easy to catch in simple describing formulae. The 
difhculty is that there is a different dependence in scale 
on the parameters e and N in the middle of the chain 
and at the ends of the chain. A rough estimate indicates 
the existence of a zone of length ^/N at the ends of the 
chain with the typical bending over of the average (yi). 
In the middle remains a zone of length order N, with 
fairly constant correlations. This splitting up in "bulk" 
and "surface" behavior which keep each other in balance 
prevents a systematic expansion in the small parameter 
e. 



VI. DISCUSSION 

Using the DMRG technique we have determined the 
properties of the RD model for moderately long chains. 
At zero and for very small driving fields we could reach 
chains of the order of -/V ~ 100 — 150 reptons; for fi- 
nite fields the lengths are restricted to some 30 reptons. 
This regime is far outside the domain where exact di- 
agonalization of the reptation matrix is possible. The 
DMRG results have the advantage, over corresponding 
Monte Carlo results, of being virtually exact as long as 
the iteration method converges. Since the DMRG proce- 
dure gives simultaneously all the lengths N smaller than 
the maximum, we can accurately determine the finite size 
effects on the asymptotic large N behavior. 

We find that the renewal time r and the diffusion co- 
efficient D are strongly affected by finite size corrections 
in the regime of these moderately large N. Here we have 
shown that the large finite size corrections, characteris- 
tic for the reptation process, manifest themselves as effec- 
tive exponents for the asymptotic behavior of the renewal 
time and the diffusion coefficient. We found that DMRG 
results reveal that, while the leading correction terms, 
as given by Doi's theory fit rather well the data for 
large N, they are not sufficient to cause a crossover be- 
havior and higher order corrections need to be included. 
These finite size effects offer an explanation for the dis- 
crepancies between measurement and standard theory. 
This point can be further tested experimentally by play- 
ing with the structure of the end reptons (e.g. branching) 
and the coordination number of the embedding lattice. 
We have lumped these aspects in a "dimensionality" pa- 
rameter d, which indeed has surprising effects on the fi- 
nite size behavior. In particular we expect that chains 
with short branching ends will be mapped onto small d 



regimes, where r and D, according to our results, will 
scale with effective exponents deviating from the values 
measured so far. Experimental tests of this prediction 
will possibly provide new insight for the understanding 
of the dynamics of entangled polymer melts and concen- 
trated solutions. 

We have also established the onset of the so-called 
band inversion. For fixed driving field and increasing 
N we observe that the drift velocity goes through a min- 
imum. This is another intriguing effect contained in the 
RD model. The band inversion has been ascribed to the 
fact that long polymers in a driving field get oriented and 
that due to this fixed orientation the drift velocity rather 
increases with N than decreases as in the non-oriented 
regime. To see the ultimate asymptotic velocity longer 
chains than presently possible should investigated. 

The DMRG calculations also yield a host of detailed 
information about the structure of the reptating polymer 
as for instance the local correlation functions. We have 
plotted the average values (yi) and (yiyt+i). We find 
that the reptating chain develops a rich and delicate pat- 
tern of shapes and correlations which are not so easy to 
catch in simple describing formulae. The difficulty is that 
there is a different dependence in scale on the parameters 
e and N in the middle of the chain and at the ends of the 
chain. A rough estimate indicates the existence of a zone 
of length ViV at the ends of the chain with the typical 
bending over of the average (yi). In the middle remains 
a zone of length order N, with fairly constant correla- 
tions. This splitting up in "bulk" and "surface" behavior 
which keep each other in balance prevents a systematic 
expansion in the small parameter e. 

We are grateful to G. T. Barkema for helpful discus- 
sions and to T. Lodge for drawing our attention to the 
experimental results of Refs. ^ and 
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